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We present convergent gravitational waveforms extracted from three-dimensional, numerical sim¬ 
ulations in the wave zone and with causally disconnected boundaries. These waveforms last for 
multiple periods and are very accurate, showing a peak error to peak amplitude ratio of 2% or bet¬ 
ter. Our approach includes defining the Weyl scalar +4 in terms of a three-plus-one decomposition 
of the Einstein equations; applying, for the first time, a novel algorithm due to Misner for computing 
spherical harmonic components of our wave data; and using fixed mesh refinement to focus resolu¬ 
tion on non-linear sources while simultaneously resolving the wave zone and maintaining a causally 
disconnected computational boundary. We apply our techniques to a (linear) Teukolsky wave, and 
then to an equal mass, head-on collision of two black holes. We argue both for the quality of our 
results and for the value of these problems as standard test cases for wave extraction techniques. 
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I. INTRODUCTION 

At the present time, several ground-based gravitational 
wave detectors, using laser interferometry, are already 
or very near to operating Qil , and the NASA-ESA 
space-based antenna LISA is scheduled to be launched 
in the next decade 3- These experiments should pro¬ 
vide the first direct probe of strong-field gravitational 
physics. The data analysis needs of these observations, 
however, require accurate waveform templates for use in 
matched-filtering algorithms. While the early and late 
stages of a merger process can be treated analytically 
using post-Newtonian and perturbation theory, respec¬ 
tively, the highly dynamical merger period can only be 
understood with the full, non-linear Einstein equations. 
In this latter regime, numerical relativity is essential 0. 

Computing waveforms from three-dimensional numer¬ 
ical simulations has, however, proven challenging SHE, 
j§j]. One of the primary problems is that there are a va¬ 
riety of length scales in a typical problem with sources. 
For a binary black hole collision, for example, the dy¬ 
namics of the merger scale with the masses and spins of 
the black holes, while the waves generated by those dy¬ 
namics have length scales that are one or two orders of 
magnitude larger. In addition, the boundary conditions 
applied at the edges of the computational domain tend to 
generate spurious ingoing radiation (which may or may 
not satisfy the full Einstein equations), which can easily 
contaminate wave signals unless the boundary is causally 
disconnected from the place and time at which the signal 
is desired. In addition, one must choose a stable form of 
the 3+1 Einstein equations, supplemented by an appro¬ 
priate set of gauge conditions, in order to build a code 
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(cf. Ref. Q). Only recently has it become possible to mix 
all of these components into a single simulation in three 
dimensions. (See Refs. 00 for very recent examples.) 

In this paper, we bring together a variety of tech¬ 
niques and technologies to successfully simulate dynam¬ 
ical spacetimes; to compute, via the Weyl scalar 'F 4 , a 
gauge invariant measure of gravitational radiation; to an¬ 
alyze that radiation quantity through spherical harmonic 
decomposition; and to demonstrate that our approach 
not only converges very well, but that it also is very ac¬ 
curate. For our sample linear problem, a Teukolsky wave 
propagating on a flat background, we see peak error to 
peak amplitude ratios on the order of 0.4%, and for our 
sample non-linear problem, the head-on collision of two 
equal mass black holes, the same ratio is on the order 
of 2% or less. These waveforms maintain this level of 
accuracy for several periods. 

While we present our results, we also wish to lay out 
some well defined examples for future use as testbeds 
for wave extraction schemes, much in the spirit of the 
“Apples with Apples” tests for numerical relativity 0- 1 
To this end, we attempt to abstract away, where possi¬ 
ble, from the details of our Einstein solver, and focus on 
( 1 ) the higher level details of how we define and com¬ 
pute gravitational radiation quantities given numerically 
evolved solutions to the Einstein equations, (2) the de¬ 
scription of a logical progression of test cases for vali¬ 
dating results, and (3) the enumeration of analytic so¬ 
lutions, symmetries, and cross-checks between test cases 
that lead, when taken together, to a high level of confi- 


1 While this paper was in final preparation, Ref. was released 

as a preprint, and, indeed, used very similiar test cases for wave 
extraction. Ref. Il111 focuses primarily, but not exclusively, on 
extractions via the Zerilli-Moncrief formalism rather than the 
Newman-Penrose formalism that we use here. 
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dence in both the methods used and the results obtained. 

To these ends, we carefully describe our particular ap¬ 
proach to the problems of computing and analyzing grav¬ 
itational waveforms in Section QH In that section, we 
describe our method for computing the Weyl scalar H /4 
from our 3+1 numerical data, describe our particular al¬ 
gorithm for computing spherical harmonic components of 
the Weyl scalar, and briefly describe our Einstein solver. 
In Section mu we bring this machinery to bear on the 
Teukolsky wave spacetime. This is a linear problem with 
an analytic solution, which we use to validate our code 
and our methods. We then turn to a non-linear problem, 
the head-on collision of a binary black hole system in 
Section ITvl which more closely resembles astrophysically 
interesting sources and tests our methods on a non-linear 
problem. We conclude with some discussion of our re¬ 
sults in Section El The appendices contain more details 
about the spherical harmonic decomposition algorithm 
and about the Teukolsky spacetime. 

II. METHODOLOGY 

Computing gravitational waves from numerical simu¬ 
lations of the Einstein equations requires combining a 
variety of nearly independent mathematical formalisms 
and numerical techniques into a single code. The meth¬ 
ods that we present, implement, and test here are not 
specific to any particular formalism or to our code. They 
are made explicit in Section HEU where we define grav¬ 
itational radiation via the Weyl scalars in the Newman- 
Penrose formalism, and in Section iHBl where we dis¬ 
cuss the need and a method for computing spherical har¬ 
monic components of radiation data. In Section urn 
we briefly discuss our particular Einstein solver, with 
which we carry out numerical simulations to validate our 
method. 


A. Weyl Scalars 

We use the Newman-Penrose formalism to compute 
gravitational radiation quantities. In this formalism, one 
chooses a tetrad of four null basis vectors, which are con¬ 
ventionally labeled Z a , n a , m a , and m a . (Note that, al¬ 
though our code uses a 3+1 decomposition of the Einstein 
equations, indices on Newman-Penrose quantities are al¬ 
ways understood to run over four dimensional spacetime 
indices a = 0,1, 2,3.) Of these vectors, Z° and n a are real 
and point along “outgoing” and “ingoing” directions, re¬ 
spectively. The vectors m a and fh a span “angular” direc¬ 
tions and are complex conjugates of each other. These 
vectors are always chosen to satisfy the orthogonality 
conditions 

l a m a = n a m a = 0 . ( 1 ) 

In addition, we impose the normalization conditions 
l a n a = —1 and m a fh a = + 1 , which are used by most, 


but not all, authors (cf. Ref. T3| l. 

Given a tetrad, tensor quantities can be recast as sets 
of coordinate scalars by projecting tensor components 
onto the basis. Applying this procedure to the Weyl ten¬ 
sor * 2 Cp qrs yields five complex scalars. This is particularly 
useful for gravitational wave studies since one of the five, 

\I /4 = —C pqrs n p fh q n r fh s , (2) 

represents outgoing gravitational waves. Indeed, the pri¬ 
mary goals of this paper are to construct 4*4 in some sam¬ 
ple, numerically evolved spacetimes, and to demonstrate 
that such computations are accurate characterizations of 
the wave content of the spacetime. 

Since, in our code, we have only 3+1 quantities on a 
single time slice from which to compute the four dimen¬ 
sional Riemann tensor, we follow Ref. 0 and write 

^ Rijkl = Rijkl + 2Ki[kKl]j (3a) 

{i) R 0j ki = -2(K j[kJ] +T™ [k K l]m ) (3b) 

^Rojoi = Rji ~ KjmK™ + KKji (3c) 

which express four dimensional quantities on the left 
hand sides in terms of 3+1 quantities on the right . 3 Here 
the indices are spatial ( i,j,k,l,m = 1,2,3), K i3 is the 
extrinsic curvature tensor for a spatial slice as embedded 
in the full spacetime manifold, and I\ is its trace. The 
symbols T”),, Rijkl, and Rjl stand for, respectively, the 
connection coefficients, Riemann tensor, and Ricci tensor 
associated with the three-dimensional spatial metric on 
the slice. 


B. Spherical harmonic decomposition 

In our simulations, we would like to be able to ex¬ 
tract the spherical harmonic components of gravitational 
waves. We find this very valuable when analyzing the 
data because 

1. Numerical noise tends to have higher angular fre¬ 
quency than genuine wave signals, and is therefore 
filtered by the decomposition process. 

2. A priori knowledge about symmetries in the data 
or dominant modes associated with physical pro¬ 
cesses allow important checks on the plausibility 
of numerical solutions (especially when exact solu¬ 
tions are not available). 


- In vacuum, which is the only case that we consider here, the 

Weyl tensor is equal to the Ri em ann tensor. 

3 There are two errors in Ref. Hi associated with what we call 
Eq. @ First there is a factor of two difference between Eq. J^bJ, 
which is correct, and the corresponding equation in Ref. Il4j . 
Second, Eq. 0 in this paper properly reflects the fact that the 
left hand sides refer to the four dimensional Riemann tensor. 
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3. Some characteristics of gravitational radiation, 
such as quasi-normal modes for instance, are best 
understood in terms of spherical harmonic compo¬ 
nents. 

In practice, however, we face a technical problem in 
computing spherical harmonic components of our data. 
While our data is stored on a cubic grid, the spherical 
harmonic components 

r) = j s Yim(0 , r, 6 , <j>)dD. (4) 

of a general (complex valued) function $ in a spherical 
harmonic basis with spin-weight s are defined by an in¬ 
tegral over a sphere of radius r. Computing the integral 
in Eq. 0 requires some type of interpolation from grid 
points to points on that sphere. 

One could solve this problem by the straightforward 
method of first interpolating the grid function $ to points 
on the sphere, and then using some numerical approxi¬ 
mation to the integral in Eq. 0 . This process of inter¬ 
polation followed by integration would then need to be 
performed at each time, at each radius, and for every 
function for which the spherical harmonic components 
are desired. 

We adopt a different algorithm due to Misner. Fol¬ 
lowing Ref. [lil? we smear the surface integral in Eq. 0 
into a volume integral over a shell of half-thickness A, 
and create an orthonormal basis for functions on this 
shell by combining the spherical harmonics in the angu¬ 
lar directions with the first N Legendre polynomials P n 
in the radial direction. This approach, which we describe 
in more detail in Appendix El has the advantage that, 
given the parameters A and N, one need only compute a 
relatively small number of weights to carry out the vol¬ 
ume integral, and that these weights depend only on the 
grid structure and the radius of extraction. This means 
that the weights can be computed once at the beginning 
of a simulation and stored. Further interpolations are not 
needed, and the weights are valid for all functions. (The 
weights do depend on the grid structure, so they would 
needed to be recomputed after grid structure changes in 
a simulation using adaptive mesh refinement.) 

The question of how to choose the input parameters 
A and N was addressed in Refs. El E3|- Under the 
assumption (motivated by the analysis in the references) 
that A is chosen proportional to the grid spacing, the 
parameter N controls the order of convergence of the 
decomposition algorithm, while the parameter A controls 
the size of the error at that order. This analysis leads to 
the 

Rule of Thumb: Choose N just large 
enough to ensure that the error term propor¬ 
tional to A is not of leading order in the grid 
spacing h. Choose A just large enough to 
safely resolve Pn on the shell. 

For second order accurate codes such as ours, Refs. El 
E2 suggest choosing N = 2 and A = 3/i/4 based on this 


rule and empirical experiments. In cases where the shell 
passes through multiple refinement regions, the grid spac¬ 
ing of the coarsest grid through which the shell passes 
should be chosen to ensure that Pn is resolved on both 
sides of the interface. We follow these suggestions for all 
work presented here. 


C. The Hahndol Code 

Our code is a fully three-dimensional, non-linear Ein¬ 
stein solver. We use a fairly standard formulation of 
the 3+1 Einstein evolution equations known as BSSN 
(T1 flg| ; our particular implementation was described in 
detail in Ref. [§(|. Because the formulation of the equa¬ 
tions is not of primary interest here, and because the 
basics of the BSSN system are widely known, we do not 
focus on these equations. 

The code uses second-order accurate finite differenc¬ 
ing to approximate spatial derivatives and the iterative 
Crank-Nicholson method M to integrate the evolution 
equations forward in time. In both our previous and cur¬ 
rent work, we have verified second order convergence in 
our results. 

Our parallel code uses the PARAMESH libraries (2^ to 
handle domain decomposition and inter-processor syn¬ 
chronization. In addition, these libraries enable us to 
use non-uniform grids to focus computational resources 
in specific areas of the computational domain. Although 
the libraries support, and we look forward to using, the 
ability to adaptively modify the grid structure during the 
course of a simulation (adaptive mesh refinement), we 
currently fix our grid structures in advance using a priori 
estimates of where to focus resolution (fixed mesh refine¬ 
ment). Our code shows 92% of optimal scaling, up to 256 
processors, when simultaneously doubling the number of 
processors used for a simulation and keeping the number 
of data points per processor constant |2(j . 

Previous studies showed that the Hahndol code is able 
to propagate linearized gravitational waves (defined in 
terms of metric components using Teukolsky’s solution 
fUfl) across mesh refinement boundaries |24j, and that 
the same code can handle strong, dynamical potentials 
moving across refinement regions |20|. In this paper 
we again focus on wave propagation, this time using 
a more formal definition of gravitational radiation, the 
Weyl scalars in the Newman-Penrose formalism, and on 
analyzing such data in a meaningful way. 


III. TEUKOLSKY WAVES 

Our first problem studies Teukolsky’s solution [§5) 
to the linearized Einstein equations, which represents a 
weak gravitational wave propagating through space. The 
linear nature of the initial data makes this an excellent 
first test problem for two reasons. First, the there is a 
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analytic solution for all times, allowing a direct calcu¬ 
lation of numerical errors for convergence and accuracy 
studies. Second, because this solution consists of a lin¬ 
ear wave in the initial data, we are able to extract our 
waveforms at small radii and short evolution times. This 
second fact is extremely useful for debugging algorithms, 
and also allows higher resolutions for a fixed problem size 
since the volume of the computational domain need not 
be as large as it would in problems (like the head-on col¬ 
lision described in Section 113 in which the waves are 
generated by non-linear sources and must propagate to 
the wave zone before being extracted. 


A. Analytic Preliminaries 

Although the Teukolsky solution is well-known, we 
summarize it here for completeness and to establish no¬ 
tation. The general form of the spacetime metric 

ds 2 = — dt 2 + (1 + Af rr )dr 2 

+ 2 Bf r grdrdO + 2Bf r( j,r sin 9drd<p 

+ (l + Cfgg + r 2 d9 2 (5) 

+ 2 (A — 2 C)fg ! j > r 2 sin ddOdcj) 

+ (l+ Cf$+Af$)r 2 sin 2 9dcj> 2 

is given in terms of angular functions fij and functions 
A , B, C. The angular functions for the l = 2, m = 0 
(spin-weight —2) case that we consider here are given by 
Eqs. (IB II) in Appendix^] The remaining functions, given 
by Eqs. (IB2II . are written in terms of a free generating 
function F. 

We follow Choi et al. |3 i n choosing 

F{x) = ^ e -* 2 /A 2 (6) 

as the exact form of the generating function, where the 
free parameters A and A represent the amplitude and 
the width of the wave respectively. The natural length 
unit in the problem is A, and we consistently choose A = 
2 x 1CC 6 A. Moreover we take an equal superposition of 
an ingoing wave and an outgoing wave, both centered at 
the origin, for the initial data. This particular choice has 
a moment of time symmetry that allows us to set the 
extrinsic curvature tensor to zero in the initial data. 

Choosing F of the form in Eq. © gives a waveform 
with oscillations but of essentially compact support. This 
is ideal for testing codes with boundaries (both inter¬ 
nal and external) because it makes clear when the wave 
passes through those boundaries, and allows one to eas¬ 
ily detect any reflections that occur due to poor interface 
conditions; cf. Ref. p4$ . 

For our 3+1 code, we use this analytic solution to gen¬ 
erate the initial data, and evolve forward in time using 
the gauge conditions lapse a = 1 and shift j3 l = 0. 


We use the Kinnersley tetrad [2^ 


l 

n 

m 


1 


(r 2 + a 2 , A, 0, a) 


2£ 


(r 2 + a 2 , - A, 0, a) 


1 


V2(r + ia cos 9) 


(iasin0,O, l,*csc0) 


(7a) 

(7b) 

(7c) 


to extract the waves in this linear wave problem. 4 
Eqs. Q give the tetrad in terms of the Boyer-Lindquist 
coordinates for a Kerr black hole. In these equations, 
A = r 2 — 2 Mr + a 2 , Y, = r 2 + a 2 cos 2 6 , M is the mass 
of the Kerr black hole, and a is its spin. Because Boyer- 
Lindquist coordinates are compatible with the coordi¬ 
nate system used for the Teukolsky metric (Eq. ©) when 
M = a = 0, we find this coordinate expression for the 
tetrad completely acceptable for this problem. 

For this relatively simple spacetime, one can compute 
analytically the value of the Newman-Penrose quantity 


vk 4 


sin 2 6 
16 


d 2 C d 2 A 
— 12 „ _ + 6 „ „ + r 


dt 2 


dt 2 



d 3 A 

~dt 3 


which represents outgoing gravitational radiation, 
ing that 


( 8 ) 

Not- 


- 2 * 20 (0)0) = \j7^ sin 2 9 (9a) 

= (y O o(0, 0) - -^20(0,0)) (9b) 

it is clear that, as claimed, Eq. m is a pure l = 2, m = 0 
mode. 

Note also that Eq. (l9EI) provides a way to compute 
a spin-weight —2 spherical harmonic component from 
the usual (spin-weight 0) spherical harmonic components. 
We make use of this in the results presented here because 
our current implementation of the Misner algorithm com¬ 
putes only spin-weight 0 components, which suffices for 
our current work. In the near future we will have the ca¬ 
pability to compute spin-weight —2 components directly. 


B. Numerical Results 

We evolved initial data specified by the Teukolsky so¬ 
lution described in Section [ill Al on a domain with outer 
boundaries at 32A using a mesh with fixed refinement 
levels having a “box-in-box” structure centered on the 
origin. The boundaries were at 2A, 4A, 8A, and 16A; 
neighboring regions differed in resolution by a factor of 


4 The first component of m a should have the sin© term in the 

numerator as in Eq. ITcl : the corresponding equation in Ref. Q 

is incorrect. 
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two, with the finest regions surrounded by coarser re¬ 
gions. We need to apply a matching condition at the in¬ 
terfaces between grids of different resolutions; we use the 
algorithm described in Ref. 0. We exploit the symme¬ 
try of the solution to reduce our computational burden 
by evolving only one octant of our data and using ap¬ 
propriate symmetry boundary conditions to mimic a full 
grid. In order to compute convergence factors, we ran the 
code at innermost resolutions A/24 and A/48. The wave 
propagates outward from the origin, crossing the refine¬ 
ment boundaries in turn. We then attempted to compute 
the spherical harmonic components of that wave at five 
distinct spherical “detectors” located at r = 3A, 4A, 5A, 
6 A, and 7A. 

To help visualize how the refinement boundaries and 
the extraction radii are related geometrically, we find 
it useful to draw two dimensional extraction maps , as 
shown in Fig. |U The extraction maps make it clear that 
the shells used to compute the spherical harmonic com¬ 
ponents generically pass through multiple refinement re¬ 
gions, especially since in three dimensions the corners of 
the cubic refinement regions tend to pass through the 
spherical extraction shells. The extraction maps shown 
in Fig. n correspond to the innermost (r = 3A) extraction 
radius and the r = 5 A extraction radius in two sample 
resolutions. (Extraction maps for nearly all extraction 
radii and two resolutions appear in Ref. |l0|.) Note that 
only the innermost refinement regions are shown in the 
maps. In each case, the entire map was surrounded by 
additional (coarser) refinement regions. 

We present our numerical waveforms in Panel A of 
Fig. El which shows the l = 2 , m = 0 component of 
T 4 for a single wave as computed at our five distinct 
extraction radii. The amplitude scales approximately as 
1 /r, which is the correct behavior to leading order in 1 /r, 
and the shape of the wave remains consistent as it moves 
outward. The consistency between the waveforms as ex¬ 
tracted at the various radii is demonstrated convincingly 
in Panel B, in which we have scaled-up the waveforms 
by a factor or r and shifted them in time to t = 3A. 
Assuming a unit speed of propagation, this should, up 
to sub-leading terms in 1 /r, result in overlaying wave¬ 
forms. Indeed the agreement demonstrated in Panel B is 
striking. 

In evaluating the effects of the refinement boundaries 
on the spherical harmonic decomposition algorithm, one 
should bear in mind that no two of the extraction radii 
have the same geometric relationship with the refinement 
boundaries (cf. Fig.|^), and yet they nonetheless generate 
perfectly consistent results. 

Since we know the analytic solution for this test prob¬ 
lem, we were able to make a detailed convergence study 
with two resolution levels. We verified that not only 
the waveforms themselves but also the extracted spheri¬ 
cal harmonic components converge at second order, even 
with the extraction spheres passing through the refine¬ 
ment boundaries, and that the solution is highly accu¬ 
rate. Fig. El shows the errors, scaled by appropriate pow- 


Extraction Map: (24, 3, 0.25) 



Extraction Map: (48, 5, 0.25) 
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FIG. 1: Selected extraction maps. Each map shows the com¬ 
putational grid in a coordinate plane and is labeled with 
a triple of numbers (N,R, A), which indicate the number 
of (cell-centered) grid points across one coordinate direction 
in the coarsest region, the extraction radius, and the half¬ 
thickness of the shell, respectively. Lengths are measured 
in units of A. Note that the shells generically pass through 
multiple refinement regions, especially since, in three dimen¬ 
sions, the corners of the cubic refinement regions tend to 
poke through the spherical extraction shells. The three cir¬ 
cles drawn on each graph show the extraction sphere (red) 
and the edges of the finite-thickness shell (blue) around the 
sphere used in the Misner algorithm. 


ers of two to indicate second-order convergence, of the 
numerical results as compared to the analytic solution. 
It is clear from the figure that we do have second order 
convergence. Note from Fig.Elthat our peak error in our 
highest resolution simulation is approximately 6 x IIP' 
for a waveform with peak amplitude of approximately 
1.5 x 10 -4 (cf. Panel A of Fig. Elk giving a fractional 
error of approximately 0.4%. 
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FIG. 3: A convergence plot for the Teukolsky wave simulation. 
The lines are the errors in the l = 2, m = 0 (spin-weight 
—2) component of the simulations at the two resolutions as 
compared to the analytic solution. 


use the excision algorithm of Ref. [3Qi] to avoid evolving 
the portions of our grid interior to the (initially separate) 
apparent horizons of our initial data. We neither allow 
our excision regions to move nor switch to a single, larger 
excision region after a common horizon forms. Although 
one or both of these techniques will likely be desirable in 
future simulations, we find the simple approach sufficient 
to extract convergent waveforms in this case. 

Once we generate our initial data, we evolve it in our 
3+1 code using the 1+log evolution equation 


FIG. 2: The l = 2, m = 0 (spin-weight —2) component of 
the Teukolsky wave, as computed at five different radii for 
the highest resolution run. The waveform is preserved up to 
the leading order 1/r scaling. Panel A shows the raw data. 
Panel B shows the same data, scaled-up by r ex t and shifted 
in time to t = 3A, the location of the inner-most extraction 
radius. 


IV. BLACK HOLE HEAD-ON COLLISION 


- = —2aK (10) 

for the lapse a, and the hyperbolic Gamma-driver evolu¬ 
tion equations 


9/3* 

~dt 

dB l 

~dt 


ui 
~W 


(lla) 

(lib) 


We also tested our techniques on a well-studied non¬ 
linear problem, th e eq ual-mass, head-on collision of two 
black holes [T?l O, HfJ |2ij. This problem is more closely 
related to the astrophysical problems of interest to us 
and, consistent with our goal of testing our techniques, 
allows a variety of symmetries that we can use as checks 
on our numerical solutions. 


A. Analytic Preliminaries 

In our head-on collisions, we place two black holes, 
each of mass M/2, on the z axis at a coordinate distance 
of 1.1515M from the origin. The initial data is generated 
by the puncture prescription of Ref. [28j, which is a gen¬ 
eralization of Brill-Lindquist initial data J29|. We then 


for the shift /3*. In these equations, K is the trace of the 
extrinsic curvature tensor for our slice, T* is the confor¬ 
mal connection function of the BSSN evolution system 
(cf. Ref. |Ti^l. and B l is defined by Eq. (Ill al to make 
the gauge evolution equations first order in time. The 
quantity 

«r) = l + (12) 

i= 1 1 1 

is the Brill-Lindquist factor used to generate puncture 
initial data for N black holes with masses Mi and posi¬ 
tions r). (In our case, N = 2 and M\ = M 2 = M/2.) 
In our work, we choose the parameter 77 = 2.8/M , and 
set a = 1 and /3* = 0 at the initial time. These gauge 
conditions were first studied in Ref. Q- 
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For the head-on collision we do not use the Kinner- 
sley tetrad, since we have no way of knowing that the 
coordinate expressions in Eq. 0 are appropriate to this 
numerically evolved spacetime. (Indeed, for a general 
problem, one would assume that they are not appropri¬ 
ate.) Instead, we construct an orthonormal tetrad from 
the numerically evolved spacetime using a Gram-Schmidt 
procedure as described in Ref. Q. 5 

Since we do not have an analytic solution for the Weyl 
scalars for the head-on case, we cannot directly compute 
the error in our solution. The symmetries of the prob¬ 
lem nonetheless provide us with several analytic checks 
on our numerical results. Before discussing these sym¬ 
metries specifically, however, it is worth making explicit 
a related point: There are three independent symmetry 
axes in this (and any axisymmetic) problem. The first is 
the symmetry axis of the physical problem, i.e. the axis 
along which the two black holes collide. The second is 
the axis with respect to which the tetrad is computed. 
(Careful examination of the Kinnersley tetrad defined by 
Eq. 0, for example, reveals this explicitly in coordinate 
expressions.) The third is the axis with respect to which 
we compute the spherical harmonics. Only if we align 
all of these axes of symmetry (conventionally along the 
2 axis) will all of our symmetry checks be true. 6 

Assuming that, as in the problem described here, all 
three symmetry axes are aligned along the 2-axis, the 
numerical solution should have the following properties: 

• Rej'Ic} is axially symmetric and is symmetric un¬ 
der the transformation 2 —> —z up to the round-off 
level. 

• Im{T4} converges to zero. 

• Viewed in a spin-weight —2 basis, the dominant 
contribution to Rej’I'.i} comes from the l = 2, m = 
0 mode. 

• Viewed in a spin-weight 0 basis, truncation level 
errors in Im{4 , 4} appear in only in modes with odd 
values of l > 3 and even values of to ^ 0. 


5 We used the Teukolsky wave as a check on our implementation 
of the Gram-Schmidt construction. In the limit of perturbed 
flatspace, the Gram-Schmidt procedure recovers the Kinnersley 
tetrad up to a trivial rescaling of l a and n a . Comparisons of 
waveforms extracted from Teukolsky wave spacetimes with the 
two choices of tetrad match very closely when this rescaling is 
taken into account. 

6 In addition, because, for example, the inner product 
f — 2 ^ 22 ( 0 , 0)oV2(0) 4>)dQ yf 0 for all l > 2, it is impossible 
to construct the spin-weight —2 components from spin-weight 
0 components in problems lacking axial symmetry, unless one 
computes spherical harmonic components at all values of l in the 
spin-weight 0 basis. 


B. Numerical Results 

For the head-on collision simulations, we use a box-in¬ 
box mesh refinement scheme similar to that employed in 
Section IT TT Bl We place our outer boundary at coordinate 
distance 128M, and place refinement boundaries at 4M, 
8M, 16M, 32M, and 64M (six levels total). We find 
that in the non-zero shift case, we need to modify our 
matching condition at the refinement boundaries; this 
will be the topic of a forthcoming paper. We evolve only 
one octant of our spacetime, using appropriate symmetry 
boundary conditions to mimic a full grid. Because of the 
octant symmetry conditions used, only one of the black 
holes appears in our evolved numerical grid. (The other 
is accounted for by the symmetry boundary conditions.) 
We ran simulations at three resolutions with innermost 
resolutions of M/16, M/24, and M/32 in order to per¬ 
form convergence tests. 

We excise a cubical region centered on each puncture. 
The cube has sides of length 0.23125A7. The size, shape, 
and location of the excision region remains fixed during 
the course of the run. 

In these simulations, we extract our waveforms at r = 
20M, 30M, 40M, and 50M. Since this is a non-linear 
problem, we are forced to use spheres at larger radii than 
in the Teukolsky wave simulations in order to ensure that 
the signals are being extracted in the wave zone. 

Using the data from these three resolutions, we were 
able to evaluate the convergence behavior of our wave¬ 
form. Fig. 0] shows a three-point convergence plot for 
the l = 2 , to = 0 (spin-weight —2) component of 
^4 extracted from our simulations. In the three point 
convergence graphs, we show (4 , 4 Mcd ' ) — \|/^ Hlsh )) and 
7(4 , 4 Low1 — 4 , 4 Med ^)/20, which should coincide, up to the 
effect of higher order error terms, for our second order 
accurate code, and our choice of relative grid spacings . 1 
The agreement is impressive. Panel A shows the conver¬ 
gence of the l = 2 , to = 0 component of ^4 as a function 
of time and as extracted at radius r = 20 M. Panel B 
shows the same comparison for the extraction radius at 
r = 40M. We observe that there is a slight degrada¬ 
tion of the signal at larger radii, but that this is to be 
expected because the larger radii are located in coarser 
regions of our grid. 

In Fig. 0 we show the l = 2, m = 0 components of 
^4 as extracted at all four radii. As with the Teukolsky 
wave, we show, in Panel A, the raw waveforms plotted 
as a function of time, and, in Panel B, the scaled and 
shifted waveforms. Again note the excellent agreement 


' The unusual convergence factor 

]_ = (1/3) 2 - (1/4) 2 
20 (1/2) 2 - (1/3) 2 

is consistent with our simulations, which have resolutions in the 
ratio 1/2 : 1/3 : 1/4. 









t (M) 

(B) 

FIG. 4: A convergence plot of the l = 2, m = 0 (spin-weight 
— 2 ) component of 'I /4 for the equal mass black hole head-on 
collision problem. Panel A shows the convergence at the r = 
20M extraction radius, and Panel B shows the convergence 
at the r = 40 M extraction radius. Note that, consistent with 
the coarser resolution at larger radii, there is a slight decrease 
in the agreement between the scaled errors in Panel B, but 
that the agreement is still excellent. 



t(M) 


(A) 



FIG. 5: A comparison between the l = 2, m = 0 components 
of as computed at distinct radii. The waveforms have been 
scaled-up by r and shifted to the innermost extraction radius, 
r = 20 M. In order to account for the fact that our lapse 
is not geodesic, we scale the waveforms by an approximate 
Schwarzschild radius R = r(l + M/(2r)) 2 rather than the 
simulation’s radial coordinate r, and shift in time by the light- 
travel time in a Schwarzschild background. 


between waveforms extracted at different radii manifest 
in Panel B. 

It is important to note that our scaling and shifting 
procedure is more complicated in the head-on collision 
case than it was in the Teukolsky wave case. In the 
head-on case, we need to account for the fact that that 
the coordinate speed of light is not spatially constant due 
to the presence of singularities. We estimate the appro¬ 
priate scaling factor and time shift by (1) assuming that 
we extract our waves far enough from the binary system 
that the metric is approximately Schwarzschild, (2) using 
the corresponding Schwarzschild radius 


R = r 



(13) 


as the scale factor rather that the simulation’s coordinate 


radius r, and (3) computing the time shift 

AT = [R' + 2Mln(R'/2M - l)] R , =Ro (14) 

in terms of our reference (Schwarzschild) radius Rq. The 
expression in square brackets is often called the “tortoise 
coordinate.” In the case of Fig. 0 Rq = R(20M) = 
21.0125 M. 

Eq. (TT1 was derived in Ref. (20] (among other places), 
and is simply the coordinate conversion between the 
isotropic radial coordinate r and the Schwarzschild radius 
R for a single, spinless black hole. It is, strictly speaking, 
only applicable to our simulation on the initial time slice. 
At later times, our slices differ from the Schwarzschild 
slices due to our choice of lapse a; far enough from the 
black holes, the slicings should be consistent and Eq. m 
should be a better and better approximation. 

Eq. 1141) gives the Schwarzschild coordinate time AT 
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required for an outgoing, radial, null ray to travel from 
radial coordinate i? 0 to radial coordinate R. It is also a 
good approximation to our coordinate time difference At 
in the limit of large radius. 

Noting that both Eq. m and Eq. Q are large-radius 
approximations to the appropriate scale factor and time 
shift, it is interesting to see that, although the waveforms 
extracted at all radii shown in Fig. agree very well, 
the extractions at the two farthest radii, r = 40 M and 
r = 50M, are extremely consistent. This is especially 
apparent in the first peak at t — AT ss 25 M and first 
trough at t — AT « 30 M. 

In addition to checking the convergence of our wave¬ 
forms and verifying consistency between the waveforms 
computed at different radii as shown above, we also veri¬ 
fied that the four “sanity” checks enumerated at the end 
of Section EH are satisfied, namely that both the real 
and imaginary parts of tI/4 have all applicable symme¬ 
tries, and that the imaginary part converges to zero with 
increasing resolution. The combination of all of these 
checks gives us confidence in our results. 

Since we do not have an analytic solution for this prob¬ 
lem, we are not able to directly compute the error in our 
simulation. As a first estimate of the accuracy that we 
achieve, we compute the ratio between the peak error, as 
computed by the difference between the high and medium 
resolution waveforms, and the peak amplitude of the high 
resolution waveform. In the worst case, our extraction at 
r = 50M, this ratio is approximately 2%. 


V. DISCUSSION 

We have presented a complete approach to comput¬ 
ing gravitational radiation via the Newman-Penrose for¬ 
malism given a numerically evolved solution to the 3+1 
Einstein equations. We validated our approach on the 
linear Teukolsky wave problem, for which we have an an¬ 
alytic solution. This allowed us to very carefully study 
the convergence behavior of our results, and to precisely 
measure the accuracy of the code. 

Starting from the Teukolsky wave results, and com¬ 
bined with the technologies of mesh refinement and a 
spherical harmonic decomposition algorithm compatible 
with a non-uniform grid, we have also demonstrated that 
we are able to compute highly convergent waveforms gen¬ 
erated by genuinely non-linear sources at sufficiently far 
distances to be considered in the wave zone. At the same 
time we locate the boundary of the computational do¬ 
main at a sufficient distance to causally disconnect it 
from the extraction region and still compute several cy¬ 
cles of the waveform. 

The high quality of our waveforms, taken together with 
the consistency between waveforms extracted at various 
radii and with different geometric relationships to our 
mesh refinement boundaries, also validates our particu¬ 
lar choices for algorithms. It specifically indicates that 
our treatment of the mesh refinement boundary condi¬ 


tions does not introduce significant numerical errors into 
our simulation, and that our spherical harmonic decom¬ 
position algorithm, given good data on a non-uniform 
grid, does not introduce any significant errors related to 
the relative positions of the extraction radii and the re¬ 
finement boundaries in the grid. 

The work described in this paper sets the stage for us 
to study more interesting astrophysical cases. We are 
currently extending these studies to non-equal mass col¬ 
lisions and inspirals, in which we will extend our existing 
mechanisms to compute, based on the spherical harmonic 
components of our waveforms, the energy and momen¬ 
tum radiated by these systems. We also hope that these 
results and these well-defined test cases can serve as a 
basis for future code validations and comparisons much 
in the spirit of Ref. El- 
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APPENDIX A: SUMMARY OF MISNER’S 
METHOD 

This appendix provides a more detailed look at the 
Misner algorithm for computing spherical harmonic com¬ 
ponents of a function represented on a cubic grid. Ad¬ 
ditional details can be found in Misner’s original paper, 
Ref. [l5} . A detailed discussion of the truncation error as 
a function of the algorithm’s parameters can be found in 

Refs. [HQ-l. 

Recall that the problem is to compute the spherical 
harmonic components defined by Eq. 0 , of a func¬ 
tion $ that is known only on vertices of a cubic lattice. 
(In this appendix, we suppress the spin-weight index s.) 
In order to begin, two definitions are need. First define 
a family of radial functions 

R n (r;R, A) = (“a~”) ( A1 ) 

in terms of the usual Legendre polynomials P n . Here R 
and A are parameters that will be associated with the 
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radius at which the spherical harmonic decomposition is 
desired and half of the thickness of a shell centered on 
that radius. From this, define 


Making use of this new metric, and with some further 
analysis, the approximation for the spherical harmonic 
coefficients 


Y nlm (r, 0, </,) = R n (r)Y lm (0, </>) (A2) 

which form a complete, orthonormal set with respect to 
the inner product 

(/Iff) = [ f(x)g{x)d 3 x (A3) 

Js 

on the shell S = {(r, 0, cj>) \ r € [i? — A, R + A]}. Note 
also that, because the functions R n form a complete set 
in the radial direction, 

R) = j p{r;R, A)Y) m (0, </>)4>(r, 0,(/>)d 3 x (A4) 

with 

OO 

p(r; R 1 A) = J2 Rn(R ; R, A )R n (r-, F, A), (A5) 

n—0 


R) = ^ R lrn{x; R)w x $(t, x) (A10) 

xer 

follows with 

N 

Rim(r ; R) = J2 Rn(R)Y nlm {r , 0,0) (All) 

n —0 

in terms of Y A = G ba Yb , not Ya- 

Note that one need only store the combination w x Ri m 
at points where w x ^ 0 in order to compute the spheri¬ 
cal harmonic components. This buries all the details of 
the computation in a relatively small table of numbers, 
and, since these numbers are not time dependent, this 
calculation need be done only once per simulation. 

APPENDIX B: TEUKOLSKY WAVE SOLUTION 


and that 


p(r-,R,A)=r 2 5(R — r) (A6) 

is a delta function. (Compare Eq. (IA4I) to Eq. QJ.) 

On a finite grid T, the inner product Eq. (TO l will have 
the form 


(/Is) = f{x)g{x)w x (A7) 

rrer 


where each point has some weight w x . This weight was 
given the form 


0 \r - R\ > A + h/2 

h 3 \r - R\ < A - h/2 

(A + h/2 — \R — r\)h 2 otherwise 


(A8) 

by Misner, where h is the grid spacing. Only cases with 
A > h/2 are considered. This means, roughly, that 
points entirely within the shell S are weighted by their 
finite volume on the numerical grid, points entirely out¬ 
side of the shell S have zero weight, and points near the 
boundary are weighted according to the fraction of their 
volume inside S. 

With the numerical inner product Eq. JUl), and letting 
capital Roman letters A = ( nlm) represent index groups, 
the Ya are no longer orthonormal. Their inner product 


{Ya\Yb) = Gab = Gba (A9) 


The general form of the spacetime metric for the 
Teukolsky wave solution j23| is given by Eq. JSJ. The 


angular function are 

f rr = 2— 3 sin 2 8 (Bla) 

f r e = — 3 sin 0 cos 0 (Bib) 

U = 0 (Blc) 

= 3 sin 2 9 (Bid) 

/^ = "I (Ble) 

fH = 0 ( Blf ) 

/$ = -fee (Big) 

4J = 3 sin 2 0-1 (Blh) 


for the l = 2, m = 0 case. The remaining functions 
' ir(2) 3F (i) 


A = 3 


B = — 


( F 


C = 


1 /FW 
4 \ r 


1 

3^(2) 

^3 

2^(3) 


3F 

ry 5 

6 F(B 

r 4 

9F( 2 ) 


6 F 


21F(B 
4- a -1" 


(B2a) 

(B2b) 


are written in terms of a free generating function F = 
F(t — r), which we choose to have the form given by 
Eq. ©. The notation 


forms a metric for functions on the shell. (Although a 
priori this matrix appears to be complex valued, it is ac¬ 
tually real-symmetric and sparse, cf. Refs, [la, 13- Bor 
now it suffices to follow Misner in denoting it as generi- 
cally Hermitian.) The inverse to this metric G AB can be 
used to raise indices on functions defined on the sphere. 


p( n ) 


d n F{x)~ 


dx n 


x=t—r 


(B3) 


denotes various derivatives. Taking F as a function of t — 
r corresponds to outgoing waves. To generate an ingoing 
solution, change the argument of F to t + r, and change, 
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in Eq. (IB2I) , the sign in front of all of the terms with odd contribute to the sum. 

numbers of derivatives. (Note that the description of how 

ingoing waves are constructed in Ref. [3] contains an 

error. The description here, which matches the original 

reference, Ref. [ 23 , is correct.) 

The Weyl scalar 4 /4 for this spacetime is computed 
from the definition Eq. © using the Kinnersley tetrad, 

Eq. Q, and noting that, of the twelve non-zero compo¬ 
nents of the Riemann tensor associated with the metric 
Eq. ©, only 


Rtete 


Rtere 

Rnprc/) 

Rr&r6 

Rtc/>r(j) 


3 2 . 2 d 2 C 

~2 T S1 " V 


1 ,a 2 A 

—r~- 

2 dt 2 


-r 2 sin 4 ( 


d 2 C d 2 A 


dt 2 

1 2 • 2 n d2A 

+ 2 


dt 2 


--r 3 sin 2 8 ( 3 


d 3 B 

~W 


d 3 A\ 
~dt 3 ~) 


- sin 2 8R t ete 

. 2 Q Rt<t>t<t> 

sin 8 

- sin 2 8Rte r e 


(B4a) 


(B4b) 

(B4c) 

(B4d) 

(B4e) 

(B4f) 
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